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Abstract. We provide a consistent statistical-mechanical treatment for describing 
the thermodynamics and the structure of fluids embedded in the hyperbolic plane. 
In particular, we derive a generalization of the virial equation relating the bulk 
thermodynamic pressure to the pair correlation function and we develop the 
appropriate setting for extending the integral-equation approach of liquid-state theory 
in order to describe the fluid structure. We apply the formalism and study the influence 
of negative space curvature on two types of systems that have been recently considered: 
Coulombic systems, such as the one- and two-component plasma models, and fluids 
interacting through short-range pair potentials, such as the hard-disk and the Lennard- 
Jones models. 
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1. Introduction 

When a system is embedded in a space of nonzero curvature, its long-distance properties, 
such as the thermodynamic quantities and the critical behavior, are strongly modified. 
A uniform positive curvature, characteristic of spherical geometry, has the drawback of 
leading to spaces of finite extent. Genuine long-distance properties on the other hand 
are accessible in hyperbolic spaces characterized by a constant negative curvature. The 
infiuence of the negative curvature of hyperbolic space on critical phenomena has been 
investigated both for hard spins on hyperbohc lattices[Il[2l[3lll[3[6l[71[8l[9l|10l[IIlll2l[l3] 
and in a continuum field-theoretical description [HI El flSl fT6]. 

In brief, the presence of an intrinsic characteristic length associated with the 
metric makes the critical behavior of Ising-type systems mean-field-like [H [21 HI [151 El [H] 
and drastically alters, if not suppresses, critical points in models with a continuous 
symmetry |14l [T6l HI [11] In the last decade, the thermodynamics of simple fiuids in 
the hyperbolic plane has also received special attention. This is the case for systems 
interacting through the Coulomb potential, such as the one- and two-component plasma 
models[ni [TBI [l9] and for the hard-disk fluid[20l [2ll [22]. 

In addition to being a mathematical curiosity and providing simple systems for 
studying the infiuence of space curvature on long-distance properties, fiuids in the 
hyperbolic plane have been considered in the context of jamming phenomena and 
glass formation in liquids. The rationale behind this comes from the concept of 
"geometric frustration" that describes an incompatibility between the preferred local 
order and the tiling of the whole space [231 [21] • Geometric frustration has emerged 
in the theoretical description of glasses and amorphous solids from the consideration 
of local icosahedral order in metallic glasses. No such frustration however operates 
in assemblies of monodisperse spherical particles in ordinary ("fiat") two-dimensional 
space, and Nelson and coworkers [25l [261 l2l] have proposed to curve two-dimensional 
space in order to introduce frustration and thereby mimic glasses in three-dimensional 
space. Two recent investigations have been carried out along these lines. Modes and 
Kamien[20l [21] have focused on isostatic packings in hyperbolic space, the isostatic 
property having been argued to play an important role in the physics of nearly jammed 
random packings of particles in Euclidean three-dimensional space [27j. On the other 
hand, building on a frustration-based approach of the glass transition [28l [29] . we have 
studied by molecular dynamics simulation the slowing down of the dynamics with 
decreasing temperature in a monatomic liquid: a truncated Lennard- Jones model in 
the hyperbohc plane [30]. 

In the present work, we provide a consistent statistical-mechanical treatment for 
describing the thermodynamics and the structure of fiuids embedded in the hyperbolic 
plane. The paper is organized as follows. In section [21 we introduce the statistical- 
mechanical framework for systems in the hyperbolic plane. By using a variant of the 
Green-Bogoliubov method, we relate the bulk thermodynamic pressure to the pair 
distribution function. In section [3l we extend the integral-equation approach, which 
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is commonly used in liquid-state theory, to calculate the pair correlation function in 
hyperbolic geometry. We discuss thermodynamic consistency and gas-liquid critical 
behavior. In section IU we focus on Coulombic systems, such as the one- and two- 
component plasma models. We apply the formalism developed in the previous sections 
to two cases for which the pair correlation functions are known exactly: the limit of 
high temperature and a special value of the temperature. In section [5l we present the 
numerical solution of standard approximate integral equations of liquid-state physics 
(namely, the Percus-Yevick and hypernetted chain equations) that have been obtained 
after making use of the harmonic analysis in the hyperbolic plane. We consider the 
hard-disk fluid and the (truncated) Lennard- Jones liquid and compare with existing 
simulation data. Finally, section [6] contains our concluding remarks, and some technical 
details are provided in two appendices. 



2. The generalized virial equation of state 

For an atomic liquid at equilibrium in the grand-canonical ensemble (fixed temperature 
T, total area A and chemical potential /i), the partition function is given by 

N=l ' i=l j>i 

where u{r) is the pair potential with r the appropriate interatomic geodesic distance (for 
simplicity, we only consider pairwise additive interactions), P = l/^ksT), and Zi(T,A) 
is the partition function of a single particle, i.e. the ideal gas contribution. The above 
expression is valid for any two-dimensional space with dA the appropriate infinitesimal 
surface area. It therefore applies as well to the hyperbolic plane, with however some 
subtelties that are discussed below. 

The hyperbolic plane is a two-dimensional manifold of constant negative 
curvature, whose metric in polar coodinates (r, (p) is given by 

d.^ = dr^+f^^^Vd0^ (2) 



The Gaussian curvature of is — k^. The curvature parameter k > therefore 
measures the deviation from "flat" (Euclidean) space; by an abuse of language, we shall 
sometimes refer to as the "radius of curvature". (Note that is not embeddable 
in three-dimensional Euclidean space and is usually described through representations 
such as the Poincare disk model that is a conformal projection of onto a disk of 
radius k,~^.) In polar coordinates, the differential area is expressed as 

sinh(Kr) , , , . , 

dA = ^ — -drdd). (3) 

The hyperbolic plane is of infinite extent but the exponential character of the metric at 
large distance (see Eqs. ([2]) and ([3D) induces the peculiar property that the boundary 
of a finite region of grows as fast as the total area of this region when the latter 
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increases. A simple illustration is provided by a disk of radius r: from Eq. ([3]), one finds 
that its area is equal to A{r) = 27r(cosh(/tr) — 1)//?^ whereas its perimeter is given by 
P(r) = 27rsinh(Kr)/K, so that for large r, both A{r) and P(r) grow as exp(Kr). 

As a consequence of the above property, boundary effects are never negligible, even 
in the thermodynamic limit. Generally speaking, the statistical mechanics of fluids on 
involves multiple integral over space, as in Eq. ([1]). Integration is performed on both 
the "bulk" and the "boundary region" of the system, the latter being taken as a region 
of finite width (or a width whose ratio to the linear size of the system goes to zero when 
the latter goes to infinity) near the boundary in an otherwise very large sample. The 
boundary contribution of course depends on the boundary condition imposed on the 
system, but it is never negligible compared to the bulk contribution. This is already 
true for the ideal- gas limit. The canonical partition function Zi is proportional to the 
total area A; however, the proportionality constant is not a pure bulk property and 
depends on the type of boundary condition [21]. 

In this work, we are only interested in the bulk behavior of fiuids in H^, which 
is what can be directly compared to fiuids in Euclidean space. One way to eliminate 
unwanted boundary effects in is to consider periodic boundary conditions[31J. This 
has been done for instance by Modes and Kamien for calculating the properties of a gas of 
hard disks in the low-density limit near ideal behavior pi]. We shall take here a simpler 
route that has already been followed in previous studies of spin and field-theoretical 
models in hyperbolic geometry [H HI [15]: we implicitly restrict all integrals appearing 
in any statistical-mechanical expression to the "bulk" of the sample, i.e. the fiuid far 
enough from the boundary. All resulting thermodynamic and structural quantities will 
therefore be dubbed "bulk" ones. We stress again that the present problem associated 
with the thermodynamic limit in hyperpolic geometry is a general one which is not 
related to the presence of long-range correlations between particles. (This point will be 
addressed further down.) 

From the (bulk) grand-canonical partition function in Eq. ([1]), one obtains the 
(bulk) thermodynamic pressure as 



For short-range interactions, the low-density expansion of the equation of state can 
be put in the usual form [21] . 



where p = N/A and Bj is the jth virial coefficient. When restricting the calculation to 
the bulk of the system, the second virial coefficient B2 can be expressed as 



where dA is given by Eq. ([3]) and /(r) = (e — 1) is the Mayer function. (In 
the calculation of B2 for short-range interactions, it is easy to see that one obtains the 




(5) 




(6) 
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same result by restricting the system to the bulk region and by using periodic boundary 
conditions [21].) As a result, 

B,^^ r,r'J^f(r), (7) 
Jo 

which for hard disks of diameter a simplifies to [21] 

(cosh(/ta) - 1) 

-D2 = TT — . (8j 

The above expression correctly reduces to the Euclidean result B2 = when k — >■ 0. 

To go beyond the low-density regime and describe the liquid phase, a natural 
strategy is to relate the pressure to the pair distribution function[32]. However, the 
hyperbolic geometry prevents a direct derivation through the Clausius function |32j. We 
have instead used a variant of the Green-Bogoliubov met hod [33]. To calculate the 
derivative in Eq. (jl]), we choose to perform an affine transformation of the elementary 
area element, 

dA' = (1 + OdA, (9) 

with ^ an infinitesimal parameter. Inserting Eq. ([3]) in Eq. Q yields the following 
transformation of the radial coordinate r to first order in ^: 

, ^(coshf/tr) — 1) 

r=r + ^- ■ w , • (10) 

K smh[Kr) 

The presence of curvature leads to a nonlinear transformation of the coordinate, contrary 
to what occurs in fiat Euclidean space P^. We then consider the infinitesimal variation 
of /n(S) generated by the affine transformation, 

<51n(5) = e(iV)-| J dA,j dA2p'^^\v,,r,)5u{r,2), (11) 

where p*^^-*(ri, r2) is the (bulk) two-particle density. By using Eq. ffTOl) . the infinitesimal 
increase of the potential is found equal to 

^ , , ^cosh(Kr) — 1 ,, , 

5u{r) = e }J . u'{r) (12) 

where u'{r) denotes the derivative of u{r) with respect to r. 

Taking advantage of the homogeneity and the isotropy of (far from the 
boundary), we then obtain the bulk thermodynamic pressure as 

^^^r.-^rir,(r)(.oM.r)-l)Ar), (13) 

P 1^ Jo 

where we have used that p'^^^(ri,r2) = p'^g{r), g{r) being the bulk radial distribution 
function[32]- This result generalizes the virial equation of state obtained in Euclidean 
space. The latter is recovered when taking the k — ^ limit in Eq. (fT3l) : 

= 1 — / drr g[r)u[r). (14) 

P 2 Jq 

At large distance, g{r) goes to 1, and Eq. (fT3|l shows that the pressure is only 
defined if (ir(cosh(Kr) — l)'u'(r) is finite, with R some irrelevant cut-off distance. This 
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imposes that the pair potential u{r) decays faster than exp(— /tr) when r — > oo. This 
requirement prevents using any interaction potential with an algebraic decay, contrary 
to the (i— dimensional Euclidean case where all potentials u{r) oc r~" with n > d lead 
to a well-defined thermodynamic limit (for the Coulombic interaction, see below). The 
above condition is of course satisfied for hard disks, but also for the one- and two- 
component plasma models because, as will be discussed later, the Coulomb potential 
decreases exponentially fast at large distance in the hyperbolic plane. 

In the case of hard disks, after introducing the auxiliary function y(r) = g{r)e^'^^^\ 
we can express the pressure as 

— — = 1 H — ^ /" dr{cosh(Ka) — l)y{r) — exp(—Pu{r)) (15) 
P Jo dr 

and using ^ exp(— /^^(r)) = S{r — a) finally leads to 

^ = 1 + ^(cosh(/€a) - l)gia+). (16) 
p 

At low density, the radial distribution function at contact 5'(cr+) goes to one, so that one 
recovers the virial expansion in Eq. ([5]) to first order with the second virial coefficient 
given by Eq. ([8]). 



3. Integral equations 

In liquid-state theory, approximate integral equations are a standard tool to describe 
the structure and the thermodynamics of a system. The approach is based on the 
(exact) Ornstein-Zernike equation that relates the radial distribution function g{r), or 
more precisely the so-called pair correlation function h{r) = g{r) — 1, to the direct 
correlation function c(r) which is obtained from the second functional derivative of the 
grand potential with respect to local density fluctuations [32], [35]. In the hyperbolic 
plane, the Ornstein-Zernike equation reads 

h{r) = c(r) + pj dA'h{r')c{t{r, r')) (17) 

where t{r, r') is the modulus of the displacement associated with an element of the 
hyperbolic translation group [36] . (In Euclidean space, t{r,r') = |r — r'|.) Again, Eq. 
f fT7|) should be considered as a bulk result, in which boundary effects have been removed 
when taking the thermodynamic limit. 

The Ornstein-Zernike equation is more conveniently expressed after applying a 
Fourier-Helgason transform (that generalizes the Fourier transform of the Euclidean 
space): 

h{k)=d{k)+pd{k)h{k). (18) 

Note that for an isotropic function h{r) as considered here, the Fourier-Helgason 
transform becomes a Mehler-Fock transform [37]. For more details, see Appendix A. 

The basis of the integral-equation approach for liquids is that the direct correlation 
function c(r) has a simpler structure and is shorter-ranged than h{r) and, as a 
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consequence, is a better starting point for approximations. Common approximations 
are the Percus-Yevick (PY) and hypernetted chain (HNC) closures in which the direct 
correlation function is taken asl32l 



c(r) = (1 + 7(r))(exp(-/3u(r)) - 1) (PY) (19) 
c(r) = exp(-/?M(r) + 7(r)) - (1 + 7(r))(HNC), (20) 

where 7(r) = h{r) — c(r). 

In Euclidean geometry and for hard spheres, the PY closure (the direct correlation 
function is then zero beyond a) is amenable to an analytical solution in spaces of odd 
dimensions |3 2 1 [38l [39l l40l 1411 142] and compares rather well with computer simulation 
results. Recently, a semi-analytical solution of the PY integral equation has also been 
obtained in two-dimensional Euclidean space [13| HI]: the first ten virial coefficients 
have been calculated and compared to Monte Carlo results. In addition, the PY 
approximation for hard discs has been numerically solved on a sphere [15]. 

Due to the approximations involved in the integral-equation approach, going from 
the structure of the fluid to its (bulk) thermodynamic properties depends on the route 
chosen. Once the pair correlation function is known, one may flnd the pressure with the 
help of the virial equation, Eq. (fT3l) . but also through the compressibility relation or 
that for the excess internal energy [32]. For instance, the compressibility relation reads 

pkBTxT= l + P j dAh{r), (21) 

where again only bulk contributions are considered. 

At this point, it is important to stress a peculiar feature of the Fourier-Helgason 
transform: contrary to what occurs in Euclidean space, J dAh{r) ^ /i(0). Indeed, by 
using the results of Appendix A, one flnds that 

27r /"^ 

h({S) = — / (ir sinh(Kr)P_i/2(cosh(/tr))/i(r) 







K 



K 



7^ — / dr smh.(Kr)h(r), (22) 







where P-i/2{x) is a Legendre function of the flrst kind (conical function). As a 
consequence, the bulk compressibility in is not given by the k = value of 
the Fourier-Helgason transform of the pair correlation function. On the other hand, 
neglecting the boundary effects, one still has 

pkeTxT = / V (2^) 

1 — p J dA c[r) 

so that the bulk pressure can be obtained by a thermodynamic integration according to 
1- - f'dp'p' [ dAc{r;p') (24) 



P PJo 

where c(r, p) denotes the direct correlation function at the density p. 
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For the hard-sphere fluid in Euchdean space, the Percus-Yevick integral equation 
provides a very accurate equation of state when the pressure is calculated with the 
following empirical rule [32] : 

P = (25) 

where Pc and Py denote the "compressibility" and "virial" pressures obtained from Eqs. 
f lM|) and ([H]), respectively. 

Note finally that according to Eq. ( 12T1) . a finite compressibility implies that 
J dA h{r) is finite, which in turn requires that h{r) decays at large distance faster than 
exp(— Kr). A gas-liquid critical point with a diverging compressibility may thus occur 
with a pair correlation function that decreases exponentially as exp(— /tr). Assume for 
instance that h{r) behaves as exp(— r/^) at large distance. One then finds, with again 
the restriction that the integration is only over the bulk region far from the boundary, 
that 



2 



o c I ^ + ^^3, (26) 

when L +cxd, where reg denotes regular, finite terms. The above expression diverges 
when ^ 1^. If ,^ is a regular function of T and p when ^ < k^^, this suggests 
that the critical behavior found near to C,c = ^{Tc, Pc) = is of mean-field type, 
with classical exponents. In particular, the bulk compressibility xt would diverge as 
1/ (1 — (^0^) ~ (T — Tc)~^, implying that the critical exponent 7 = 1. In addition, 
assuming that h{r) ~ exp(— Kr) at large distance at the critical point and using the 
known expression of the Fourier- Helgason transform of 1/ cosh(/tr)[37j, it is easy to 
derive that the Fourier-Helgason transform h{k) of the pair correlation function is finite 
at = 0, even if the compressibility diverges (which is a dramatic illustration of Eq. 
fl22|) ). and has a regular expansion in {k/nY when {k/n) — > 0. 

The expected mean-field nature of the gas-liquid critical point in has the 
very same origin as the mean-field character of the critical point of the Ising model 
in hyperbolic geometry [H HJ [15]: the divergence of the compressibihty/susceptibility 
is controlled by the exponential growth with distance of the differential surface area 
associated with the hyperbolic metric. An approximate description of the gas-liquid 
critical behavior will also be provided in section 15.21 



4. Coulombic systems 

We now consider systems of particles interacting through a Coulomb potential, as 
was done before by Jancovici and coworkers [1 9 [ [T8l [T7] . For a monodisperse system 
(particles with equal charges of the same sign), it is necessary to add a charged uniform 
background such that the global electroneutrality is satisfied, thereby allowing a proper 
thermodynamic limit. In the limit of point-like particles, this corresponds to the one- 
component plasma model (OCP). When considering a binary mixture of two species of 
oppositely charged particles (with charges ±g) in equal concentrations, it is possible to 
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obtain a well-defined thermodynamic limit without introducing a charged background: 
in the limit of point-particles, this is the two-component plasma model (TCP). In this 
latter case, the system is stable at high temperature, down to a temperature at which 
a collapse of pairs of opposite charges occurs. When the charged particles have an 
additional hard-core interaction and at low charge density, the system undergoes a 
Kosterlitz-Thouless transition at a lower temperature = 5^/(4^^) O HZ] • 

The Coulomb potential is formally defined as the solution of the Poisson equation, 
which is related to the Green's function of the following inhomogeneous Laplace 
equation: 

Av{r,(f)) = -27r6^^\r,(f)), (27) 

where A is the Laplacian operator and 5^"^^ (r, 0) the Dirac distribution in the appropriate 
manifold. In the Euclidean plane and for a function depending on the radial coordinate 
only, A is given by 

A.(r) = 1^ (r^) , (28) 



r dr \ dr ^ 

and the solution of Eq. ( j27l) is 

v{r) = -\n{r/L) (29) 

where L is an arbitrary length. Therefore in two dimensions, the Coulombic interaction 
potential diverges with the distance, even if the associated force remains a (slow) 
decreasing function of the distance. 

On the other hand, in the hyperbolic plane, the Laplacian acts on a function of the 
radial coordinate as 

1 9 / dv(r)\ 

^^W = -i:^7r^i^M'^^)^r^ • (so) 

smh(Kr) or \ or J 

The solution of the inhomogeneous Laplace equation is then given by 

t;(r) = -In (^tanh (^y)) , (31) 

so that the pair interaction between two point particles with charges qi and q2 at a 
geodesic distance r is equal to qiq2v{r) with v{r) given above. The curvature of the 
hyperbolic plane thus introduces a screening of the Coulomb interaction. Indeed, v{r) 
now decreases as exp(— Kr) at large distance beyond 

It is possible to derive analytical expressions for the thermodynamic quantities 
in two cases: the high-temperature limit, where the Debye-Hiickel approximation is 
asymptotically exact, and the particular value of the temperature T = q'^/{2kB), where 
a large amount of exact results have been obtained both in the Euclidean and in the 
hyperbolic plane. In the following we shall consider the two models introduced above, 
the OCP and the TCP. (To extend the investigation beyond the two cases considered 
here, one could use the HNC integral equation which is known to give a good description 
of the pair correlation function in flat space[32], but this is left for future work.) 
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4.I. One-component plasma (OCP) 

For the OCP, the charged background provides a uniform neutrahzing contribution. 
Assuming a perfect compensation amounts to replace g{r) by h{r) = g{r) — 1 in the 
expression of the virial equation of state. After taking the derivative of Eq. (I3T!) . the 
equation of state can be written as 

p K Jq sinh(fi;r) 

where, we recall, P is the bulk thermodynamic pressure. (For a discussion of the different 
definitions of pressure in the OCP in the hyperbolic plane, see Fantoni et a/. [18].) We 
shall show that the above equation of state reduces to known results in the appropriate 
limits. 

Note that generalized Stillinger-Lovett sum rules are satisfied by the OCP in the 
hyperbolic plane |17j. namely: 



(32) 



dAh{r) = -l, (33) 

! /i(r) ln(cosh(Kr/2) = -1. (34) 
J 

These rules express the fact that the system is a conductor (in which therefore electro- 
neutrality and screening hold). 

In the high-temperature limit, the Debye-Hiickel approximation becomes 
asymptotically exact and provides an analytical expression for the pair correlation 
function /i(r)[T9l [H]: 

h{r) = -/3g2Q^(cosh(/tr)), (35) 

where q is the charge of the particles and Qi, is a Legendre function of the second kind 
with an index i> given by 



^=-o + \/7 + 2^- (36) 



1 2^^^^' 

2 + V4 + ^^ 

One can check that the expression in Eq. fl35l) satisfies the two above sum rules. 

We can now combine the generalized virial equation of state that we have derived 
with the Debye-Hiickel expression for h{r). Inserting Eq. ( 135|) in Eq. (l32l) leads to the 
following equation of state: 

(3P _ ^ nWfp r^^QM 

P Ji X + 1 

For a fixed nonzero curvature (/t > 0), Eq. ( l36i) shows that z/ in the high- 
temperature limit (/5 — ^ 0); knowing that Qo{x) = | In (jzl ) and using that 



dx^^. (37) 
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one finally obtains that 

p 12^2 ^ ' 

This expression coincides with the dominant term in the high-temperature limit of the 
virial expansion derived by Jancovici and Tellez[T9]. 

On the other hand, when k ^ (faster than /?^^^), z/ ~ A/27r/3pg2/ k —>■ +oo and 
one has 

' Q^^jx) 1 

dx— — . 40 

1 + X zz/^ 

The pressure is then given by 

f — («' 

which corresponds to the result in the Euclidean limit [48] . 

We next consider the specific value of the temperature T = q^/{2kB), usually 
expressed in terms of the coupling parameter T = (3q^ = 2. After a mapping onto a 
non-Hermitian fermionic field theory, Hastings [15] has obtained the exact expression of 
the pair correlation function for Coulomb systems in the hyperbolic plane for this value 
of the coupling. For the OCP, this provides 

= (42) 

cosh(Kr/2)*'^/'/'^ 

After introducing the change of variable x = cosh(fi;r/2) and inserting Eq. (l42l) in Eq. 
we obtain 

T = ' — ^''^ 

Since Pq"^ = 2, the bulk thermodynamic pressure is finally given by 

PP ^ 2np+jl_ ^^^^ 
p Airp + 

In the low-density limit, this gives f3P/p = 1 — (2,71/ K,^)p + 0{p'^)^ which leads to the 
same second virial coefficient as found in Ref. fT9J. 

When K — > 0, one recovers from Eq. (jUj) the exact Euclidean limit 

^ ^ \ (45) 
p 2 

Conversely when 00, the pressure goes to the ideal-gas limit. This result seems to 
be quite general, namely: the influence of the inter-particle interactions vanishes in the 
large-curvature limit. 
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4-2. Two- component plasma (TCP) 

The TCP has been investigated in Euchdean space[50l[5ll[52],[53l[54]- For point particles, 
the plasma is stable in the high-temperature conducting phase, for a coupling parameter 
Pq"^ < 2. Adding a small hard-core repulsive potential allows one to obtain a stable phase 
at lower temperatures, up to a coupling jSq"^ = 4 where a Kostelitz-Thouless transition 
between a conducting phase and a dipolar phase takes place. As for the OCP, the model 
is solvable for = 2 in various geometries [55| \56\ 157] . 

The bulk thermodynamic pressure of the TCP in the hyperbolic plane is expressable 

as 

(3P TT(3q^p r , , , , , ,,(cosh(/€r) - 1) , , 

p K Jo smh(Kr) 

where /i++(r) and h (r) denote the pair correlation functions between particles of equal 

charges and particles of opposite charges, respectively. In the high-temperature/low- 
coupling limit, the Debye-Hiickel approximation provides an accurate description of the 
conducting phase. For the TCP, the linearized Poisson equation is similar to that of the 
OCP. The pair correlation functions are then given by 

h+Ar) = -h^-ir) = -^Q.M. (47) 

After inserting Eq. fH7|) in Eq. fHBl) . one obtains an equation of state in the Debye-Hiickel 
approximation which is equal to that of the OCP in the same approximation, Eq. (139|) . 

For = 2, Jancovici and Tellez[T2] have obtained the exact expressions of the 
pair correlation functions of the TCP (in the limit of very small hard-core interactions). 
However, the calculation of the equation of state via Eq. ( H6l) becomes much more 
involved than for the OCP without bringing much new insight. Therefore, we have not 
pursued in this direction. 



5. Short-range potentials 
5.1. The hard- disk fluid 

As mentioned in Introduction, the hard-disk fluid has been considered in the context 
of geometrical frustration. Approximate equations of state combining the low-density 
expansion and a description of close or nearly close packings have been proposed [20l [2T| 
|22]. 

We have solved the PY integral equation, which is known to be generally better 
than the HNC one for hard spherical particles [32j, for a large density range. Fig. [1] 
shows the resulting radial distribution function g{r) for different values of the curvature 
at an area fraction (p = p ^^(^°^'^(^^^""/^)~^) ~ q 55 (recall that the freezing area fraction in 
the Euclidean plane is around 0.7[SH])- We have also plotted the Euclidean counterpart 
corresponding to k = 0. It is striking that for a range of curvature parameter, na 
between and 0.5, all curves are essentially superimposable: in this range, the pair 
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Figure 1. Radial distribution function g{r) obtained from the PY equation at an 
area fraction (j) = 0.55 for various values of the curvature: na = 0, 0.15, 0.5, 1.5. Only 
for K = 1.5, namely a radius of curvature smaller than the particle diameter, the 
structure displays significant deviations from the Euclidean case. The inset zooms in 
the region of the first minimum and second maximum: g(r) for Ka = 0.15 and for the 
Euclidean case are indistinguishable, whereas for Ka = 0.5, the extrema are slightly 
less pronounced. 

distribution function g{r) is rather insensitive to curvature. An evolution is on the 
other hand clearly seen for stronger curvatures, e.g. na = 1.5. 

From the PY correlation function, we have calculated the equation of state through 
Eqs. f fMl) . fl2^ . and (!25|) . The results are illustrated in Fig. [2l Similarly to the 
Euclidean case, the compressibility route gives a higher bulk thermodynamic pressure 
than the virial route, and the linear combination given in Eq. ( l25l) lies in between. 
Comparison with the molecular dynamics data of Ref. [201 121] and with the equation 
of state proposed in Ref. [22] is performed in Fig. [3l At low area fractions simulation 
data and predictions share the same behavior controlled by the second virial coefficient. 
However, the situation deteriorates for even moderate area fractions. (Note that the 
number of atoms in the simulation cell with periodic boundary condition is very small, 
always less than 10.) 

It is interesting to investigate the local order of the hard-disk fluid. The geometric 
frustration induced by the negative curvature describes the impossibility of extending 
the local order present in the dense fluid to tile the whole space. This is true for a nonzero 
but small enough frustration for which the local order remains that of the Euclidean 
plane, i.e. the hexagonal order. As stressed by Rubinstein and Nelson [26], there is a 
series of increasing curvature parameters Ka for which heptagonal, octagonal, etc., local 
order can tile space without frustration. Two points however are worth stressing. First, 
such regular tilings are only a few of all possible tesselations of the hyperbolic plane. 
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Figure 2. PY equation of state of the hard-disk fluid in versus area fraction cj). (a) 
na = 0.15: the upper curve corresponds to the compressibility route, the lower curve 
to the virial route, and the middle curve to the empirical rule in Eq. (j25p . (b) PY 
equation of state obtained from Eq. for different curvatures: na — 0.15, 0.7, 1.06 
from bottom to top. 
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Figure 3. Equation of state of the hard-disk fluid in versus area fraction </>: 
comparison between the PY result combined with Eq. ((25|) (for kct — 1.06) and the 
simulation results of Ref. f2T] for kct = 1.060 (green squares) and na — 1.062 (blue 
diamonds). The dash-dotted line is the prediction of Ref. 22J and the dashed line 
corresponds to the equation of state truncated at the second virial coefficient. 



Modes and KamienpOt [21] have rather focused on so-called isostatic regular packings, 
i.e. tilings of in which each atom has 4 neighbors (the {4, 5} tesselation for instance 
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can be realized as a close packing of disks for a curvature intermediate between the two 
curvatures studied in their simulations: kct = 1.0612). However, such isostatic packings 
have a small density compared to the heptagonal {7,3}, octagonal {8,3}, etc., tilings 
considered in Ref. [26]. Second, the local order in the equilibrium hard-disk fluid has 
not actually been investigated as a function of curvature. 

Our study sheds some light on the latter question. From the knowledge of the 
radial distribution function, we have calculated the average coordination number Z via 
the formula 



where am corresponds to the distance between particles when g{r) attains the first 
minimum beyond the first peak. The results are plotted as a function of area fraction 
for different values of the curvature in Fig. HI As expected, the average coordination 
number increases with density and seems to saturate (we cannot access close packing 
values with the approximate PY integral equation). It also increases with the curvature 
parameter. The values at high density are in semi-quantitative agreement with the 
prediction of Rubinstein and Nelson [26] based on a fictive "ideal" random close packing. 



a prediction that interpolates between the exact values of the tesselations: Z = 6 
for the {6,3} tiling when k, = 0, Z = 7 for the {7,3} tiling when na = 1.0905.., 
Z = 8 for the {8,3} tiling when = 1.5285.., etc.... The coordination number is 
clearly overestimated for large curvatures, e.g. na = 1.06, in the PY approximation and 
we suspect that the validity of the approximate closures quite generally gets worse as 
curvature increases. 

5.2. The truncated Lennard- Jones potential 
The Lennard- Jones pair potential. 



is commonly used to model liquids in Euclidean space. However, as shown in 
section [3l its power-law nature makes it inappropriate for hyperbolic geometry as 
no thermodynamic limit would be reachable. Physically, the (cr/r)^ term arises from 
London dispersion forces and its electrostatic origin suggests that the dependence on the 
geodesic distance should be modified in the hyperbolic plane. Indeed, even the Coulomb 
potential that (logarithmically) increases with distance in two-dimensional Euclidean 
space is found to decay exponentially in the hyperbolic plane with an intrinsic screening 
length which is provided by the radius of curvature k~^. As for the (a/r)^^ term, it is 
merely a convenient way to mimic a steeply repulsive potential between atoms and its 
algebraic character has no physical foundation. 






71 
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Figure 4. Average coordination number versus area fraction <j) for different values 
of the curvature, as obtained from the PY equation. From bottom to top: na = 
0.15, 0.7, 0.9, 1.06. The dashed horizontal lines indicate the maximum values predicted 
by Rubinstein and Nelson [55]. 

To preserve a proper thermodynamic limit in H^, the simplest procedure is to 
truncate the Lennard- Jones potential, Eq. ( l50i) . beyond some cut-off distance Vc- this 
should not significantly alter the behavior of the liquid, provided r,. is less than the radius 
of curvature This cut-off procedure is anyhow most often employed in the Euclidean 
case in order to save computer time, and this will facilitate the comparison between the 
results in Euclidean and hyperbolic planes. In the Euclidean plane, Glandt[59l |60] has 
obtained the numerical solution of the PY and HNC integral equations for the Lennard- 
Jones interaction potential. In this case, the HNC equation becomes significantly more 
accurate than the PY equation as the density increases. (A comparison with simulations 
has been performed in Ref. [61j.) 

In the present work, we have numerically solved both the PY and HNC integral 
equations for a large range of temperature and density. As in the case of the hard-disk 
fluid, we observe a very small evolution of the structure with the curvature and the 
differences between the PY and HNC closures are very small too up na ~ 0.5 (curves 
not worth showing here). 

We first focus on the region of density and temperature close to the gas-liquid 
critical point. A simple iterative procedure (Picard method) is then inefficient for 
converging towards the solution of the integral equation. It is necessary to proceed 
with a more sophiscated iterative scheme. We have combined the Newton method with 
a conjugate gradient procedure [62l [63] (which avoids the operator inversion otherwise 
present in the Newton method[64]). We have adapted this method to hyperbolic 
geometry. Details are given in Appendix B. 
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Figure 5. Inverse of the isothermal compressibility versus area fraction (p for the trun- 
cated Lennard- Jones fluid in the PY approximation near but above the critical tem- 
perature: from top to bottom , T = 0.640, 0.610, 0.580, 0.560, 0.530, 0.515, 0.490, 0.472. 
The curvarture parameter is na = 0.5. 

We illustrate these results for the PY closure and a reduced curvature na = 0.5. 
Differents isotherms above the critical temperature are shown in FigJSl The inverse of 
the susceptibility reaches a minimum at each temperature and the successive minima 
approach zero as one gets closer to the critical point. Figure [6^ displays the pair 
correlation functions at a reduced temperature T = 0.515 for several area fractions; 
the log-linear plot clearly indicates that the decay is always exponential and faster than 
exp(— Kr) (dotted straight line). A more stringent test of the long-distance behavior of 
the pair correlation function h{r) is shown in Fig. Eb , where exp(fi:r)/i(r) is plotted as a 
function of the distance r. One can see the approach to a plateau for large enough r as 
one approaches the critical point. (The slowest decay corresponds to the area fraction 
at which the compressibility is maximum.) 

We have also considered the liquid at a rather high density. Recently, we have 
performed extensive molecular dynamics simulations of the monodisperse truncated 
Lennard- Jones liquid in the hyperbolic plane (with periodic boundary conditions) |30j . 
We have shown that a small curvature indeed prevents long-range hexagonal (or hexatic) 
ordering and leads to glass formation as temperature is lowered. Most results were 
obtained for a density pa^ = 0.852. Fig. [7] displays the radial distribution function 
g{r) at an area fraction = 0.669 and a high temperature T = 3.259 for a curvature 
parameter na = 0.1 (pcr^ is then equal to 0.852). At this state point, the PY solution 
provides a slightly better description of the simulation data than the HNC equation: the 
HNC equation underestimates the structure of the liquid at short distances (see the inset 
of Fig. [T]). When the temperature is lowered to T = 1.885 {i.e., still above the freezing 
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Figure 6. (a) Log-linear plot of pair correlation function h{r) — g{r) — 1 in the PY 
approximation for T — 0.472. The behavior as a function of is non monotonous 
and the slowest decay is for (j) — 0.118, which corresponds to the maximum of the 
compressibility. The dotted line is exp(— Kr). (b) Same data multiplied by exp(Kr): a 
convergence towards a constant plateau at large r is clearly visible as one approaches 
the critical point. The curvature parameter is kct = 0.5. 

temperature in the Euclidean plane which is around T* ^ 0.75 at this density [30]), the 
HNC and PY closures are now of comparable agreement with the simulation data, but 
the maximum of the first peak is overestimated. Overall, the agreement appears good 
with either the PY or the HNC prediction. 

6. Conclusion 

We have formulated a statistical-mechanical treatment for describing the thermodynam- 
ics and the structure of fluids embedded in the hyperbolic plane, focusing on the bulk 
behavior. We have applied the formalism to the one- and the two-component plasma 
models, the hard-disk fluid and the (truncated) Lennard- Jones liquid model. The neg- 
ative curvature of space provides an intrinsic length that screens the long-distance be- 
havior of the pair correlations: this is observed both for the Coulombic interactions and 
for the vicinity of the gas-liquid critical point. The formalism correctly reproduces the 
known results in the appropriate limits and, via approximate integral equations and 
exact thermodynamic relations, it leads to predictions that compare well with existing 
simulation data. 

A major motivation for studying liquids in the hyperbolic plane comes from the 
theoretical description of jamming phenomena and glass formation in terms of geometric 
frustration. Frustrated local hexagonal order in negatively curved space provides a 
simple model to study the glass transition. Among the many theories proposed for 
explaining the glass phenomenology, several attribute the slowing of relaxation with 
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Figure 7. Radial distribution function g{r) of the (truncated) Lennard- Jones liquid 
on for a curvature parameter Ka — 0.1, an area fraction — 0.669, and two 
different temperatures: left, T = 3.259; right, T — 1.885. The circles correspond to 
the simulation data, the full (red) curve to the PY equation, and the dashed (blue) 
curve to the HNC equation. The inset zooms in on the region near the first minimum 
and the following extrema. 



decreasing temperature to the presence of a very large number, actually exponential in 
the system size, of metastable states of low energy or free energy [65 | [66 | [67] . Techniques 
borrowed from spin-glass theory such as the replica formalism have been applied to 
liquid models in three dimensions to check for the existence of such a multitude of 
metastable states, and evidence has been provided within mean-field-like liquid-state 
approximations [HHl EHl [70]. Since in some sense, the hyperbolic metric induces in 
statistical systems some features of mean-field models, one may wonder whether a 
glassforming liquid in the hyperbolic plane such as the truncated Lennard- Jones model 
also possesses a complex free-energy landscape at low temperature. The integral- 
equation approach developed in the present formalism can be combined with the replica 
method to address this point. Work in this direction is in progress. 

P.V. acknowledges B. Jancovici for useful suggestions on Coulombic systems. 



Appendix A. Harmonic analysis in the hyperbolic plane. 

In it is possible to define a generalization of the usual spatial Fourier transform, the 
Fourier- Helgason transform. For a generic function / of the polar coordinates r, 0, it is 
defined as |37] 

f{a,t) = [ drrf0sinh(r)e^"*P^i/2+,i(cosh(r))/(r,<^), (A.l) 
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where is a Legendre function of the first kind (conical function). The inverse 

Fourier-Helgason transform is given by 

firA) = / rf"tanh(7rt)e-^"* 

p-%+.,(cosh(r))/(a,t). (A.2) 

For isotropic {SO [2) invariant) functions, the Fourier-Helgason transform reduces 
to a Mehler-Fock transform [37], the dependence on a disappears and for a curvature 
parameter k one obtains 

f{k) = — / (ir sinh(Kr)P_]^/24_jii (cosh(Kr))/(r). (A.3) 
^ Jo 

The inverse transform is then equal to 

/(r) = — [ dkktanh{7i-)P,,^,.k{cosh{Kr))f{k). (A.4) 

AtT J K / "t" K 

Let us show how one can recover the Euclidean limit when k — >■ 0. Starting from the 
integral representation of the Legendre functions of the first kind|37j. 

P'xiz) = ^!^^J^^^l^ du[z + v/i^cos(n)]^e-^-, (A.5) 

one obtains that 

^_i/2+»^(cosh(/tr)) = 

1 /"^"^ k 

— / (i'u[cosh(/tr) + sinh(Kr) cos(-u)]^^^^'^*« , (A.6) 
2vr Jo 

and when k goes to zero, one finds that 
2k J „ 

where Jo{r) is a Bessel function of the first kind. Finally, the Fourier-Helgason transform, 
Eq. (lA.3p . becomes 

/■oo 

~f{k) = 27r / rdrUkr)f{r), (A.8) 
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which corresponds to the standard Fourier transform in two dimensions. 
In a similar way, the inverse Helgason transform goes to 

/W = T^ r dkMkr)~f{k). (A.9) 



27r . „ 

The Fourier-Helgason transform satisfies the convolution theorem: if / and g are 
5*0(2) invariant functions, one has 

{j^gm = f{k)m. (A.io) 

where the star denotes the convolution product [7T]. 

Finally, by using an exact Mehler-Fock transform, we illustrate the result given in 
Eq.f l22l) : the Fourier-Helgason transform of a given function / when A; — > is not equal 
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to the integral of this function in the hyperbohc plane. We consider /(r) = cosh(Kr)~". 
With the change of variable w = cosh(Kr), the integral of / is obviously equal to 

r>00 -| 

dww-" = (A.ll) 

1 + a ^ ^ 

whereas the Fourier-Helgason transform for k = is given by 

w^Py^^udw = ^^^^ w!V' (A- 12) 

T{a)^/7^ 

where t = k/n. The second-order expansion of Eg. flA.12p is equal to 

(r(f-i))^2- 

4r (a) 

(4 xl> (1, f + I) - 4 VI/ (1, 1 + I) a + v]/ (1, f + I) + 16) 2" (r (f - i))' 



16 {2a-lf y^T{a) 



+0{t^) (A. 13) 

Apart from the limiting case a oo, the integral over the whole space is different from 
the Fourier Helgason transform at = 0. In addition, the function /(r) has a definite 
integral only when a > 1, whereas the Fourier-Helgason transform at A; = is defined 
for a > 1/2. We also note that the coefficient of the second order in the expansion in t 
is always negative for a > 1/2. 

Appendix B. Solving integral equations in H^i the numerical procedure. 

For solving integral equations, the basic method is an iterative procedure, known as 
the Picard method, which can be easily adapted to systems in the hyperbolic plane as 
follows: Starting from an initial value of 7o(r) = hQ^r) — Co(r), Co(r) is obtained by 
using the closure equation, Eq.( fT9l) for PY or Eq.( l20i) for HNC We next perform the 
Mehler-Fock transform of Co(r). From Co{k), and using Eq. (fT8!) . we calculate jnew{k) 
through the Ornstein-Zernike equation, namely. 

Finally, we transform back 'jnewik) to obtain a new function 71 (r). The solution is taken 
as a barycenter of 7o(r') and •jnewir) , namely 71 (r) = a7o(?") + (1 — a)'ynew{r) where a is 
a mixing parameter. The whole procedure is repeated and we consider that convergence 
is reached when I'jnir) — 7„+i(r)| < e where e ~ 10^^. 

Note that the Mehler-Fock transform (Fourier-Helgason transform of SO {2) 
invariant functions) is implemented by using the Legendre or conical functions 
P-i/2+it{f')- These functions are available in the Gnu Scientific Library and are tabulated 
for saving computer time. The number of points in real and reciprocal space was equal 
to 1600, and integration was performed with a Simpson rule. In regions of the phase 
diagram, far from the coexistence curves and the freezing transition line, the above 
procedure converges rapidly. 
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However, in the vicinity of the critical point, the Picard method exhibits extremely 
slow convergence and a Newton method is more appropriate for solving integral 
equations. In addition, the equations of the Newton method are not solved by a 
direct inversion of the associated operator (or of the large matrix obtained after 
discretization of the functions), but by a more efficient procedure, the conjugate gradient 
methods [6^ [72] : it consists in solving the equations by using an iterative method that 
avoids the inversion of a large matrix (or an operator). The derivation of the method is 
done below for the PY closure, but it can be easily modified for the HNC closure. 

After differenciating Eq. (fT8|) and (fT9|) . and letting /(r) denote the Mayer function, 
the equations involved in the Newton method are then given by 

(1 - pc{k))6^{k) - p{^{k) + 2d{k))Sc{k) = p{^{k) + c{k))c{k) - 7(A;XB.2) 
6c{r) - /(r)57(r) = -c(r) + /(r) (1 + 7(r)) (B.3) 

where 6c{r){6c{k)) and 6'y{r){6'j{k)) are solutions of the Newton equations. 

By taking the Fourier- Helgason transform of Eq. (IB.3p . one eliminates 6c and one 
obtains an equation of the form 

A5^{k) = B (B.4) 
where A is an operator such that 

A~g{k) = (1 - pc{k))m - Pirn + 2d{k)){fg){k) (B.5) 
and B is given by 

B = [-c[k)+f{k)+(fi){k))p{^{k)+c{k))+p{^{k)+m 

The adjoint operator is given by 

A^m = (1 - pc{k))m - Pifh^g) + 2f{^g)) (B.7) 

where * denotes the convolution product. 

In order to solve the linear system, Eq. (]B.4p . the solution is built on a sequence of 
functions in mutually conjugate directions. We introduce the inner product in Fourier- 
Helgason space: 

(X,y) = / rf/t-tanh f— ^ X{k)Y{k). (B.8) 



K \ K 

Starting with the initial conditions 



Ro = A^B (B.9) 

Po=Ro (B.IO) 

the iterative procedure is given as follows: 
• Calculate the coefficient a^, 

{Rk,Rk) m UN 
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• The residual and the solution become 



Rk - akA^APk 



(B.12) 



(B.13) 



• The next conjugate direction is given by the relation 



-Pfc+l — Rk+l + l^kPk 



(B.14) 



with 




(B.15) 



The iterative procedure is stopped when -R/c+i) < e (e ~ 10^^° in practice). 

We then take the inverse Fourier-Helgason transform of the solution 6'^{k), and with 
the help of Eq. (lB.3p . we derive 5c(r), and finally 5c{k) by a Fourier-Helgason transform. 

We next update the functions 7 and c in both real and reciprocal spaces : 
7i = 7o + ^1 and Ci = cq + 5c. The inner product b = {B,B) (where B is given 
by Eq.( ]B.6[) ) is calculated by using the updated functions 7 and c. So, as long as 6 > e' 
with e' ~ lO"*^, the iterative procedure is repeated. 
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